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Abstract 

The quasi-Monte Caxlo method for financial valuation and other 
integration problems has error bounds of size 0((log N)^N~'^), or even 
0{{logN)^N~^^'^), which suggests significantly better performance 
than the error size 0{iV"^/^) for standard Monte Carlo. But in high 
dimensional problems this benefit might not appear at feasible sample 
sizes. Substantial improvements from quasi-Monte Carlo integration 
have, however, been reported for problems such as the valuation of 
mortgage-backed securities, in dimensions as high as 360. We believe 
that this is due to a lower effective dimension of the integrand in those 
cases. This paper defines the effective dimension and shows in exam- 
ples how the effective dimension may be reduced by using a Brownian 
bridge representation. 

1 Introduction 

Simulation is often the only effective numerical method for the accurate val- 
uation of securities whose value depends on the whole trajectory of interest 
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rates or other variables. Standard Monte Carlo simulation using pseudo- 
random sequences can be quite slow, however, because its convergence rate 
is only OiN"^^^) for N sample paths. Quasi-Monte Carlo simulation, using 
deterministic sequences that are more uniform than random ones, holds out 
the promise of much greater accuracy, close to 0{N~'^) in optimal cases. Ran- 
domized versions of quasi-Monte Carlo simulation can in some cases bring 
the typical error close to 0(iV~^/^). 

This dramatic improvement in convergence rate has the potential for sig- 
nificant gains both in computational time and in range of application of sim- 
ulation methods for finance problems. An optimistic reading of the results 
suggests an effective squaring or even cubing of the sample size N. Large im- 
provements have in fact been found in a number of earlier studies [1, 11, 18], 
which were all motivated by the results of Paskov [17] on mortgage backed 
securities. 

Quasi-Monte Carlo simulation is not a magic bullet, however. The asymp- 
totic error magnitudes are the ones it is "close to" above, multiplied by 
log(A^)*^, where A; depends on the dimension s of the simulation. In high 
dimensions these powers of log(iV) do not become negligible at any computa- 
tionally possible sample size. This loss of effectiveness has been documented 
for a series of test problems in [6, 7, 8]. When simulations are cast as inte- 
gration problems the resulting integral is often of very high dimension (e.g. 
dimension 360 for a mortgage of length 30 years), so any loss of effectiveness 
at high dimensionality can affect them. 

Our first goal in this paper is to reconcile two apparently conflicting 
truths. The first is that quasi-Monte Carlo is not much better than Monte 
Carlo in high dimensions with practical sample sizes. The second is that 
quasi-Monte Carlo has been seen to far surpass Monte Carlo in some high 
dimensional examples. It is our view that success in high dimensional prob- 
lems tells us more about the integrand than about the method of integration. 
Some high dimensional integrands are indeed amenable to quasi-Monte Carlo 
simulation. Integrands of low "effective dimension", which we define in two 
ways below, are of this type. Our second goal is to give an example of a 
financial simulation, in which one can reduce the effective dimension of an 
integrand, thereby making quasi-Monte Carlo much more effective. 

The outline of this paper is the following: Section 2 gives a brief in- 
troduction to quasi-random sequences and their properties, including the 
Koksma-Hlawka inequality which is the basic estimate on integration error 
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for quasi-Monte Carlo. The dependence on dimension and the character of 
two-dimensional projections of quasi-random sequences is also discussed. Sec- 
tion 3 introduces some useful decompositions of integrands, and uses them to 
define two notions of the effective dimension of an integrand. The mortgage- 
backed security problem is formulated in Section 4. Our main technical tool 
for formulating the problem with reduced effective dimension is the Brow- 
nian bridge representation of a random walk, which is described in Section 
5. Computational results for the mortgage-backed security problem are pre- 
sented in Section 6. Conclusions are discussed in Section 7. 

We are grateful to Spassimir Paskov and Joseph Traub for a number of 
discussions and for providing us with their quasi-random number generator 
FINDER. 

2 Quasi- Random Sequences, Discrepancy and 
Integration Error 

2.1 Basic properties 

The origin of the improved accuracy of quasi-Monte Carlo methods is the 
improved uniformity of quasi-random sequences. Figure 1 shows two plots, 
each of 4096 points in two dimensions. The top is a pseudo-random sequence 
and the bottom is a quasi-random (SoboP) sequence. In the pseudo-random 
sequence there is clumping of points, which limits their uniformity. The 
cause of this clumping is that, since points in a pseudo-random sequence are 
(nearly) independent, they have a certain chance of landing very near to each 
other. The constructions used for points in a quasi-random sequence, on the 
other hand, prevent them from clumping together. 

The uniformity of a sequence of points in the s-dimensional unit cube 
P = [0, 1]* can be measured in terms of its discrepancy. This is defined by 
considering the number of points in rectangular subsets of the cube. For a 
set J C and a sequence of N points {xn)n=i define 

n=l 

Here Xj is the characteristic function of the set J, and m( J) is its volume. 
If E* is the set of subrectangles with one corner at (0, 0, . . . , 0), then the star 
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Figure 1: 2-D Projection of a Pseudo- Random Sequence (top) and a Sobol' 
Sequence (bottom) 



discrepancy is defined as: 

D*^ - sup (2.1) 

JeE* 

Some other discrepancies do not treat the origin differently from the other 
corners of/*, and there are counterparts to the L°° definition given above. 
See Niederreiter [10] for some of these other definitions and Hickernell [4] for 
some more recent generalizations. 

The importance of discrepancy can be seen from the Koksma-Hlawka 
inequality for integration error. For the integral of a function / on the 5- 
dimensional unit cube, the simulation based estimate of the integral is 

M/) = ^E/(^n) (2.2) 

n=l 

and the integration error is 

^Nif) = ;^ E /(x„) - /^^ fix) dx. (2.3) 

n=l 

The Koksma-Hlawka inequality says that 

\er,{f)\<V{f)D*^ (2.4) 

where V{f) is the variation of /. 

In one dimension, the variation is V(f) ~ Jq \df\. The definition in higher 
dimension is more complicated. Define for all fc < s and all sets of k integers 
1 < Zi < ^2 < . . . < ifc ^ 5 the quantity 

y('=)(/;n,...,ifc) = 
The variation of / (in the sense of Hardy and Krause) is defined as 

k^l l<ti<t2<-<lfc<5 

The Koksma-Hlawka inequality (2.4) should be compared with the for- 
mula for root-mean-square error of Monte Carlo integration using a random 
sequence. If the sequence x„ is uniformly distributed on then 

EMfn^^ - a(/)iV-i/2 (2.5) 



dti^ • • • dti^ 



dti^ • • • rfttfc 



5 



in which E is the usual expectation and cr(/) is the square root of the variance 
of / given by 

aU)=i^f^{m-I?dx)"^ (2.6) 

with f = fdx. 

The error magnitudes (2,4) and (2.5) are similar in that the bound is 
a product of one term depending on properties of the integrand function 
and a second term depending on properties of the sequence. The Koksma- 
Hlawka inequality is an absolute bound, which is more satisfying theoretically 
than (2.5), an equality in expectation which holds only probabilistically. For 
practical purposes the preference is reversed. Each factor in (2.4) is incredibly 
hard to compute, whereas the Monte Carlo variance can be estimated from 
the same data used to compute /. Furthermore the Koksma-Hlawka bound is 
an inequality that is only tight for a worst case function /, whose fluctuations 
are exquisitely matched to the discrepancies in the sequence (xn), while the 
Monte Carlo variance estimates the error for the actual / being sampled. 

The infinite sequence {xn)^=i is said to be quasi-random if 

D;,<c{logN)^N"^ (2.7) 

in which the constant c and the logarithmic exponent k may depend on the di- 
mension s. For integration by quasi-random sequences, the Koksma-Hlawka 
inequality says that the integration error is size O {{log N)^N~^), which for 
large N is much more accurate than standard Monte Carlo simulation. 

Examples of quasi-random sequences have been constructed by Halton, 
Faure, Sobol', Niederreiter and others. For a comprehensive discussion see 
the monograph of Niederreiter [10]. 

Although the variation of / requires s derivatives of /, we have found 
in practice that only a minimal amount of smoothness of / is needed for 
effectiveness of quasi-Monte Carlo integration. For problems in which / is 
discontinuous, however, the improvements of quasi-Monte Carlo integration 
are diminished. 

The effectiveness of quasi-Monte Carlo in high dimensional problems was 
studied in [7, 8]. They showed that the faster convergence rate for quasi- 
Monte Carlo is generally lost for problems of high dimension. The simplest 
evidence for this conclusion is seen in the discrepancy as a function of N 
for various values of the dimension s. For small dimensions, the discrepancy 
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appears to be 0{N~^), ignoring logarithmic factors, for all N. For large 
dimensions, the discrepancy behaves initially like 0(A^"^/^) as for a random 
sequence, taking on an apparent 0{N~^) rate only for very large values of N, 

The transition value of TV appears to be grow exponentially with the 
dimension. This has not been rigorously proved, but for quasi-random se- 
quences based on nets as described below, it is in keeping with the definition 
of nets. In high dimensions, unless one uses a very large number of points, 
quasi-random sequences are no more uniform than random sequences. Thus 
the Koksma-Hlawka bound does not imply any advantage for quasi-Monte 
Carlo with moderate values of N and large dimension s. On the other hand, 
we have found almost no problems for which quasi-Monte Carlo is worse than 
standard Monte Carlo. 

It is difficult to evaluate the uniformity of a sequence in a high dimensional 
space. A necessary but not sufficient condition for uniformity is uniformity 
of low dimensional coordinate projections of the sequence. The Sobol' se- 
quence used in the computations below (based on the sequence generated by 
FINDER) has excellent one dimensional projections, and in the 360 dimen- 
sional case we test, many, but not all, of its two dimensional projections are 
very uniform. In Section 3 we show how such partial uniformity is sufficient 
for the task of integrating functions of low effective dimension. 

The one and two dimensional projections of (xn) are easily graphed. 
Three dimensional projections can be investigated, with some difficulty, with 
dynamic graphics tools, such as XGobi [21]. The graph at the bottom of Fig- 
ure 1 shows a "good" pairing of dimensions using Sobol's 2^^ and 3''^ dimen- 
sions with his recommended starting values. A "bad" pairing of dimensions is 
presented in Figure 2, which shows two higher dimensions (following Sobol's 
convention for associating dimension with generating polynomial) based on 
the polynomials -h a:^ + a:^ + + 1 and x^ + x^ -h + + x^ + x + 1 and the 
starting values (1,3,5,11,3,3,35) and (1,1,7,5,11,59,113), respectively. (The 
polynomials correspond to dimensions 27 and 32 respectively in FINDER. 
However, the starting values are different.) Although this non-uniformity 
could go away if these starting values were changed, we have found that 
this type of non-uniformity is fairly typical of some of the two-dimensional 
projections of high dimensional Sobol' and Halton sequences. 

The bad behavior seen in the second plot of Figure 2 can be explained in 
terms of filling in holes. If 8192 (2^^) points are used, the plot looks almost 
identical to what is shown for 4096. However, the next 8192 points fall only 
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4096 Points of Sobol Sequence 




Dimension 27 

Figure 2: 2-D Projection of Sobol' Sequence 
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where the gaps appear. Thus by TV = 16, 384, the projection plot is almost 
perfectly uniform. The problem is that the cycle for filling in such holes can 
be too long. 

2.2 mj5)-nets and (f, 5)-sequences 

The quasi-random points we consider are, or are based on, (t, m, 5)-nets and 
{t, s)-sequences. The definitive reference on this topic is Chapter 4 of the 
monograph by Niederreiter [10]. 

An elementary interval in base 6 is a set of the form 

for nonnegative integers a^, kj with aj < b^^ . For f > 0, a sequence of iV = 6"* 
points Xn is a (i, m, 5)-net in base b if every elementary interval J in base 6 of 
volume has Rn{J) = 0, Given m, 5, and 6, the smaller t is, the better. 

An infinite sequence (x„)^i is a (t, 5)-sequence in base b if each finite 
sequence ixn)^Ab^i^ for ^ > 0 is a (t, m, 5)-net in base 6, for all m > 0. 

The net property becomes relevant for m > t, that is > b^~^^. Below this 
value of any sequence of points in even identical points, are consistent 
with the net property. The smallest N at which the net property constrains 
some fully 5 dimensional elementary subinterval (one with all kj > 0) is 6*"*"^. 
The asymptotic rate for the discrepancy of nets should therefore start at 
around A'' = 6*+*. 

The most widely used constructions of (i, 5)-sequences are due to Sobol', 
Faure and Niederreiter. The Sobol' sequences are (i, 5) -sequences in base 
6 = 2. For s = 360 the value of t is quite large for Sobol' sequences, in the 
thousands, according to Niederreiter (personal communication, 1996). The 
Faure sequences are (0, 5)-sequences in base 6 where 6 > 5 is a prime number. 
For 5 = 360, the smallest value of 6 is 361. 

For Faure sequences the net property starts to be relevant at iV = 361, 
but the asymptotics are not relevant until = 361^^^ > 10^^^. For Sobol' 
sequences the net property starts to be relevant at TV = 6* and the asymp- 
totics are relevant at N = b^'^\ Even if t is as small as 1000 these two values 
are larger than 10^^ and 10^°^ respectively. 
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Recent constructions of Niederreiter and Xing (1996) hold the promise 
of (t, 5)-sequences in base 2 with i = 5, a significantly smaller value than 
previously known possible. Even here though this suggests that the net 
property becomes relevant at N = 2^^^ > 10^*^^ and the asymptotics become 
relevant at TV = 2^20 > iQ2ie 

Despite the impractically large values of mentioned above, quasi-Monte 
Carlo works well on some high dimensional integrands, using realistic sample 
sizes. The reason is that low dimensional views of the quasi-Monte Carlo 
points can have small discrepancy and some integrands are dominated by 
low dimensional structure, as described in Section 3. For example, the points 
in Figure 1 appear to be a (t, m, 2)-net with a much smaller value of t than 
holds for all 360 dimensions of the points. 

2.3 Scrambled (t,m, s)-nets and (t, 5)-sequences 

A scrambled net is a hybrid of quasi-Monte Carlo and Monte Carlo methods, 
in which a sequence is randomized. These are defined and analyzed 
in a sequence of papers by Owen [14, 15, 16]. See also Hickernell [3]. The 
randomization is carefully constructed in the base b so that (t, m, 5)-nets 
map onto (t, m, 5)-nets, (t, s)-sequences onto (t, 5)-sequences, and each Xn 
individually has the uniform distribution on As a result InH) becomes 
a random variable with mean / and at least quasi-random accuracy. 

By replicating the randomization one can estimate the sampling variance 
statistically from the same data used to estimate the integral. The variance 
of /^(/) is o{N-'^) for any / with f f{x)'^dx < 00. 

The precise improvement over Monte Carlo depends on properties of /. 
If / is well approximated by sums of characteristic functions of elemen- 
tary intervals of large volume, then V{lN{f)) decreases rapidly. In par- 
ticular if / is smooth enough that h{x) — d^f{x)/dx^ is Lipschitz contin- 
uous, \h{x) - h{x')\ < B\\x - x'\f for some B > 0 and /3 e (0,1] then 
V{lN{f)) = 0{N-^ {log Ny-^) for scrambled nets on TV points. This pro- 
vides for errors of 0(A/'~^/^(log7V)(*~^^/^) in probability, a better rate than 
that achieved by unscrambled nets. The extra AT'^/^ factor of accuracy may 
be attributed to cancelation of some integration errors, which does not hap- 
pen with deterministic sequences. 

The asymptotic rate can be expected to commence a,t N = 6*"^*, as for 
unscrambled nets, though with realistic sample sizes, significant benefits can 
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be seen on integrands of lower effective dimension. 

3 Integral Decompositions and Effective Di- 
mension 

Suppose that we can write 

m = E /m(^)- (3-1) 

m=0 

Then / = E^Lo/m and /^(/) - E^=o^iv(/m), and e^(/) = E^lf^o eiv(/m) 
in obvious notation. 

In decompositions like equation (3.1) we will customarily arrange that 
= / is constant in x and that = 0 for m > 1. Then 

M 

eNif) = E eiv(/m) (3.2) 

m=l 

because /iv(/o) = ^- If / fm{^)fk{^)dx = 0 for m ^ fc, the decomposition 
(3.1) is said to be orthogonal. 

In simple Monte Carlo simulation, the Xn are taken independently from 
the uniform distribution on [0, 1]^. For an orthogonal decomposition, cr^(/) = 
Em=i^^(/m)- Many Monte Carlo methods have the effect of changing the 
sampling variance of / to 

1 M 

E ^m^Vm) (3.3) 

m=l 

for some "gain" constants > 0. Getting every Vm < I implies an im- 
provement over simple Monte Carlo. Even when the are not all smaller 
than one, a method in which is small whenever cr^ is large can be very 
effective. 

A case in point is antithetic sampling. Write / — /o + /i + where 
/o(x) = /, h{x) = (fix) + /(I - x))/2 - f and f2ix) = (/(x) - /(I - x))/2. 
Here l—xis interpreted componentwise. This is an orthogonal decomposition 
into constant, even and odd (symmetric eind antisymmetric) parts. Antithetic 
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sampling involves taking the points Xn and = 1 — a;„ in pairs, for 1 < n < 
N/2. One can show that in antithetic sampling r2 = 0 and Fi = 2. Antithetic 
sampling can be anywhere from half as good as simple Monte Carlo (when 
/2 = 0) to infinitely better than simple Monte Carlo (when /i = 0). When it 
is thought that erf » <7i, as for nearly linear functions, antithetic sampling 
becomes attractive. 

In an analysis of variance (ANOVA) decomposition, M = 2* — 1, and there 
is one term in (3.1) for each subset of the s components of x. The empty set 
corresponds to /q. The term for each subset is a function only of the com- 
ponents within the subset. It is convenient to replace the labels 0, 1, . . . , M 
by subsets u C {1, 2, . . . , s}. Thus f{x) = E« fu{x)^ The ANOVA decompo- 
sition is an orthogonal one, so a^{f) = ^u<^u where = a^(/u). See Owen 
[12] for definitions, Takemura [22] for some history of this decomposition and 
Hickernell [4] for some recent generalizations. 

Let gt{x) = Z)|u|=t /wW) for 0 < ^ < 5. Then gt describes the part of / 
that is exactly t dimensional and J2d=i 9d describes the part that is at most 
t dimensional. The variance of gt is cr'^igt) = J2\u\=t^i- example, 
a^{gi) = 0.99cr^(/), then we can conclude that 99% of the variance of / is 
due to the components of X taken one at a time. In such a case we might say 
that / is eflFectively one dimensional. Likewise, if a^(5i) + a^(p2) + ^^^ios) is 
close enough to we might consider / to be effectively 3 dimensional. For 
such an /, a set of points with good uniformity in every triple of variables will 
produce an accurate integral estimate, even if the points are not particularly 
uniform in some quadruples of variables. 

Definition 3.1 The effective dimension of in the superposition sense, is 
the smallest integer ds such that Z)o<|u|<ds ^^ifu) ^ 0.99a^(/). 

This notion of eflFective dimension diflFers from the one implicitly used in 
[17], which we state formally as: 

Definition 3.2 The effective dimension of f , in the truncation sense, is the 
smallest integer dx such that I^uc{i,2,...,dT} ^^(/w) — 0.99a^(/). 

The threshold 0.99 is an arbitrary choice, and one might well prefer 
other values in some settings. For each d = 1,...,5 we can define how 
d dimensional / is, in the senses above, by the ratios ^o<\u\<d^u/^^ 
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5Dwc{i,2,...,d} ^u/^^- Clearly the fraction of variance that is d dimensional is 
higher in the superposition sense than in the truncation sense. 

To distinguish the two definitions, consider the function / = Si=i ^^^i- 
For this function ds{f) = 1 and drif) = s. In general, arranging for all 
ds dimensional projections of ixn)n=i to ^^^tve low discrepancy ensures an 
accurate simulation. Alternatively, arranging for the first dr components 
of ixn)n=.i to have low discrepancy will suffice. The dimension ds does not 
depend on the order in which the input variables are indexed, while dr does. 

If / is nearly one dimensional in the superposition sense, then Latin hyper- 
cube sampling will be an extremely effective simulation technique. In Latin 
hypercube sampling (McKay, Conover and Beckman [5]) the j'th component 
of Xn is Xnj = (7rj(n) — Unj)/N where ttj is a (uniform) random permuta- 
tion of the integers 1 through N and the Unj are independent {7[0, 1] random 
variables. Latin hypercube sampling stratifies each individual dimension, 
but imposes no higher dimensional stratification. In fact, a Latin hypercube 
sample is a scrambled (0, 1, 5)-net in base N, Stein [20] shows that Latin hy- 
percube sampling is essentially equivalent to having Ft, = 0 for subsets u of 
cardinality \u\ = 1 and = 1 for |u| > 1. The variance contribution of each 
fu with |u| = 1 is o{N~^) and if that fu is smooth the variance contribution 
is in fact 0{N~^). We use Latin hypercube sampling below to investigate 
the fraction of the variance that comes from one dimensional parts of an 
integrand. 

Some orthogonal array sampling schemes [12, 13, 23] balance all margins 
up to order t, randomize the higher ones and have F„ = l|u|>t- These should 
be effective on integrands that are, or are nearly, of effective dimension t or 
less, in the superposition sense. 

For non-randomized quasi-Monte Carlo methods, the decomposition (3.1) 
does not lead to a variance interpretation, but we may still write 

M M 

Mf)] < E Mfu)\ < E ^N,uK.(/u), (3.4) 

N>i |ul>i 

where is the discrepancy of the |u|-dimensional points obtained by keep- 
ing only those components {Xn)n^i in u, and 14(/u) is the variation of /« taken 
as a |u|-dimensional function. We believe that many successes of quasi-Monte 
Carlo methods on high dimensional problems can be attributed to a low ef- 
fective dimension of the integrand, in one or both of the senses above. In 
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such cases, arranging for small values of D*pf^ to coincide with the few large 
values of al and large values of D^^^, if any, to coincide with small values of 
al (which ordinarily implies small V{fu)) should produce accurate results. 

As N increases, more and more of the terms in equation 3.4 should switch 
from the Monte Carlo rate iV-V2 ^ quasi-Monte Carlo rate N-'^{\ogN)^. 
While this is taking place, there may be an appearance of an asymptotic rate 
better than N~^^^. If however, the very highest dimensional term is nonzero, 
then one can expect that the errors are bounded below by a small multiple 
of A^~^/^ until N reaches the sometimes astronomical sample sizes ft*"*"^ as 
described in Section 2.2. 

Scrambled net sampling can be analyzed in terms of a square wavelet 
decomposition applied to each term in the ANOVA decomposition, [15]. The 
result is that the contribution of fu decreases very rapidly after TV = b*"^!"!. 
This explains how the scrambled net variance is o(iV~^). Large gains can be 
expected if / has low enough dimension in the superposition sense. 

Classical Monte Carlo techniques often replace the integrand / by another 
one with the same integral and a smaller variance. The ANOVA decompo- 
sition above suggests two more ways to improve integration: reducing the 
effective dimension, in either of the two sense above. In high dimensional 
problems, reducing both sorts of effective dimension pays off. It is, for exam- 
ple, easier to arrange low discrepancy for all three dimensional projections of 
16 variables than for all three dimensional projections of 360 variables. We 
believe that the Brownian bridge discretization in Section 5 succeeds for these 
reasons on the Mortgage- Backed securities problem introduced in Section 4. 

4 Mortgage-Backed Securities 

Consider a security backed by mortgages of length M months with fixed 
interest rate zq, which is the current interest rate at the beginning of the 
mortgage. The present value of the security is then 

PV = E{v) 

M 

= E(J2^kmk) (4.1) 
fc=i 
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in which E is the expectation over the random variables involved in the 
interest rate fluctuations. The variables in the problem are the following: 





= discount factor for month k 




= cash flow for month k 


ik 


= interest rate for month k 


y^k 


— fraction of remaining mortgages prepaying in month k 


Tk 


= fraction of remaining mortgages at month k 


Ck 


= (remaining annuity at month k)/c 


c 


= monthly payment 




= an A'^(0, a) random variable. 



This notation follows that of Paskov [17], except that our corresponds to 

his au-k-hi- 

Several of these variables are easily defined: 

rrik = cr kiil ~ Wk) -f WkCk) 

M-k 

Ck = + 

Following Paskov, we use models for the interest rate fluctuations and the 
prepayment rate given by 

ik = Ko e^* ik-i 

= Kj^e^^^'-'^^^io (4.2) 
tw/t = /C^i + K2 arctan {K^ik + ^4) 

in which Ki^K2^K^yK/^ are constants of the model. The constant Kq = 
g-«r2/2 jg chosen to normalize the log-normal distribution; i.e., so that E{ik) = 
io. The initial interest rate Zq is an additional constant that must be specified. 
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In this study we do not divide the cash flow of the security among a group 
of tranches, as in [17], but only consider the total cash flow. Nevertheless, 
the results should be indicative of a more general computation involving a 
number of tranches. 

The expectation PV can be written as in integral over with Gaussian 
weights 

5(0 = (27^a^)-'/'e-fV2-^ (4.3) 

This is transformed into an unweighted integral by a mapping ^ = G{x) with 
G'{x) = y(^), which takes a uniformly distributed variable x to an N{0,a) 
variable ^. The formula for PV is 



Note that in quasi-Monte Carlo evaluation of an expectation involving a 
stochastic process with M time steps, the resulting integral is M dimensional. 
The parameter values io, Ki, K2, K^, K4^,a used in [17] are proprietary, 
and not known to us. In the numerical study below, we have used two sets 
of values for them. The first set, chosen to be plausible to us, turned out, 
by inspection to be very nearly a linear function of the Gaussian increment 
random variables ^k- We refer to this set of parameters as the "Nearly Linear 
Example" . For this case, the parameters are 



In the second example, the "Nonlinear Example" , the parameters are 



(io, Ki, K2, Ks, K4, a^) = (.007, .04, .0222, -1500.0, 7.0, .0004) (4.6) 



This example also has a large linear component but the effect is less extreme. 
While perhaps less plausible, it represents a more challenging integrand. The 
monthly interest rate io corresponds to a yearly rate of 8.4%. The variance 
in interest rate increments leads to yearly fluctuations of size .5%. In the 
Nearly Linear Example, the prepayment rate is nearly linear in the interest 
rate, in the range of interest; whereas for the Nonlinear Example, the pre- 
payment rate has a step increase when the interest rate falls much below io- 
In both examples the length of the loans is taken to be 30 years (M = 360). 



PV 




(4.4) 



(io, Ku i^2, K^, a^) = (.007, .01, -.005, 10 



I .5, .0004) 



(4.5) 
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5 Brownian Bridge 



Since Brownian motion is a Markov process, it is most natural to generate 
its value b{t + A^i) as a random jump from a past value b{t) as 

b{t -f Ati) = b{t) + u (5.1) 

in which u is an N{0, 1) random variable. On the other hand, the value 
b{t -f- Ati) can also be generated from knowledge of both a past value b{t) 
and a future value b{T = t + Ati + ^^2), with 0 < AU, according to the 
Brownian bridge formula 

b{t + Ati) = ab{t) + (1 - a)b{T) + cr/ (5.2) 

in which 

a = Ati/(Ah-^At2) 

c = yJaAt2 (5.3) 

Note that aAi2 < Ati, so that the variance of the random part of the Brow- 
nian bridge formula (5.2) is less than that in (5.1). 

The standard method of generating a random walk = ab(kA) is based 
on the updating formula (5.1). The initial value is yo = 0. Each subsequent 
value yfc+i is generated from the previous value using formula (5.1), with 
independent normal variables 

Another method, which we refer to as the Brownian bridge discretization 
can be based on (5.2). Suppose we wish to determine the path 2/0, • • * ? VMj 
and for convenience assume that M is a power of 2. The initial value is j/o = 0. 
The next value generated is t/M = o\/MAt v^. Then the value at the mid 
point 2/M/2 is determined from the Brownian bridge formula (5.2). Subsequent 

values are found at the successive mid-points; i.e. yuj^^ Vzm/a^ Vm/s^ The 

procedure is easily generalized to general values of M. 

Although the total variance in this representation is the same as in the 
standard discretization, much more of the variance is contained in the first 
few steps of the Brownian bridge discretization due to the reduction in vari- 
ance in the Brownian bridge formula. This reduces the effective dimension 
of the random walk simulation, which increases the accuracy of quasi-Monte 
Carlo. Moskowitz and Caflisch [9] applied this method to the evaluation 
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of Feynman-Kac integrals and showed the error to be substantially reduced 
when the number of time steps, which is equal to the dimension of the cor- 
responding integral, is large. Since the mortgage-backed securities problem 
described above depends on a random walk, and can be written as a dis- 
cretization of a Feynman-Kac integral, we were naturally led to apply the 
Brownian bridge discretization to this problem, 

6 Numerical Results 
6.1 Nearly Linear Example 

The value PV for this example was calculated to be 131.78706. The mean 
length of a mortgage in this case is 100,9 months and the median length is 
93 months. The Monte Carlo variance of PV is 41.84 and the variance in 
antithetic sampling is 0.014. This suggests that the function is very nearly an 
odd function of the Gaussian increments. In fact, solving 41.84 — al-^al with 
0.014 = 2(j?, provides the rough estimates al = 0.007 and cr| = 41.833 so 
that the odd or antisymmetric part of this integrand provides about 99.98% 
of the variation. 

Similarly the variance in Latin hypercube sampling is about 0.0155 from 
which we find that roughly (41.84 - .0155)/41.84 = 99.96% of the variation 
comes from one dimensional structure. This function is effectively one di- 
mensional in the superposition sense, and it is nearly antisymmetric. The 
percentages quoted above are based on ratios of sampling variances and may 
not be exact, but both of these findings agree with what we found by nu- 
merical inspection: this function is very nearly linear in the Gaussian incre- 
ment variables. Application of Latin hypercube sampling to the antithetic 
integrand leads to only a slight decrease in variance compared to the non- 
antithetic case, which we interpret to mean that the multidimensional part 
of the integrand is not predominantly odd. 

We now describe the accuracy of various integration methods for this 
problem as a function of A^, the number of paths. For each of these results, we 
present the root-mean-square of the error among 25 computations. For Monte 
Carlo methods the 25 computations are statistically independent, at least to 
the level possible using pseudo-random numbers. The Sobol' calculations for 
each of the 25 runs and for each value of N was computed using different, 
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non-overlapping subsequences of the Sobol' sequence. Because they are not 
a sample from any population, the root mean square error presented is the 
difference between the values obtained and a "gold standard" obtained from 
using the quadratic terms of a Taylor series expansion of the integrand about 
the origin in the Gaussian coordinates as a control variate and N 3.2 
million. The results are plotted in terms of error versus N, both in log base 
10. 

In these plots, for the antithetic computations, N refers to the number 
of times the antithetic integrand {f{x) 4- /(I - x))/2 is evaluated. This 
corresponds to 2N function evaluations. Plotting the antithetic runs versus 
2N would be more appropriate when function evaluation is the dominant 
cost and plotting versus when generating the Xn is the dominant cost. 
We don't attempt to plot cpu time on the axis as this can depend on how 
efficiently a method is implemented. 

First, we perform straightforward Monte Carlo evaluation, with results 
plotted in Figure 3. The top curve shows results from Monte Carlo using 
standard pseudo-random points, with the error decreasing at the expected 
rate of iV~^/^. The second curve shows a dramatic improvement using the 
360 dimensional Sobol' sequence (generated with part of the code FINDER). 
In a separate calculation (not plotted), it was found that if only the first 
50 dimensions (using the standard discretization) were taken to be quasi- 
random and the rest pseudo-random, the size of the error decreased slightly 
compared with the purely random case, and the apparent convergence rate 
remained 7V~^/^. So, in the truncation sense, the dimensionality is not below 
50. 

The third curve in Figure 3 shows the results of combining the quasi- 
random sequence with antithetic sampling. This leads to a sizeable reduc- 
tion in the error size as the dominant anti-symmetric part has been removed. 
However, the improved convergence rate, characteristic of low-dimensional 
quasi-Monte Carlo methods, also disappears. Finally, reference lines for Latin 
hypercube sampling and antithetic random Monte Carlo are shown. These 
both effectively remove the one dimensional linear elements of the integrand, 
with antithetic variates kilUng it off exactly, while this error decreases like 
0{N~~^^^) for the Latin hypercube case, leaving the remaining errors to con- 
verge at the standard N~^/^ rate. The quasi-random error, while outperform- 
ing simple random Monte Carlo on the one dimensional elements, still can 
only achieve the 0{N~'^) rate in the optimal case. Both Latin hypercube and 
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antithetic random outperform the antithetic quasi-random sequence slightly. 
This may be because the quasi-random sequence is performing at worse than 
the Monte Carlo level on the higher dimensional part of the integrand. 

A scrambled (0, 360)-sequence in base 361 with N a small multiple of 
361 behaves like Latin hypercube sampling. When N nears 361^ = 130321 
some bivariate sources of variance disappear. For intermediate multiples of 
361, some of the bivariate effects are reduced by a multiple On 
this function, the results for the scrambled (0, 360)-sequence were essentially 
that same as for Latin hypercube sampling, and thus they are not plotted in 
Figure 3. When applied to the antithetic integrand, the scrambled sequence 
showed a slight improvement over the antithetic Latin hypercube sampling; 
however, no significant gains are achieved in either case with antithetic sam- 
pling due to the increased computation time required. 

Our results are consistent with the results of Paskov [17, 18] and Ninomiya 
and Tezuka [11]. They contradict the observation that the effectiveness of 
quasi-Monte Carlo is lost in high dimensions. However, we argue below that 
the improvement is almost entirely due to improved integration of the one 
dimensional parts of the integrand. 

Next we consider the Brownian bridge version of the integrand. The 
reformulated integrand has the same mean and the same variance as the 
original, but more of the structure is packed into the first few dimensions. 
This should help the Sobol' sequence because it would only need to have 
small discrepancy among the first few dimensions. This encoding decreases 
the variance in Latin hypercube sampling from 0.0155 to 0.00963, suggesting 
that the BB encoding has made the integrand even more inherently one 
dimensional. 

Figure 4 shows the results from Sobol' sequence integration in the Brown- 
ian bridge representation, with and without antithetic sampling. Also shown 
are reference lines for simple Monte Carlo, which is not affected by the 
change of representation, and for Latin hypercube sampling with the Brow- 
nian bridge. Again antithetic sampling does not substantially improve Latin 
hypercube sampling here. For the Brownian bridge Sobol' sequence without 
antithetic sampling, the results are essentially the same as for the standard 
representation. This is because for the dominant linear one dimensional el- 
ements, the Brownian bridge representation simply rearranges the weights 
on the elements, but the sum remains constant. Because the errors associ- 
ated with each one dimensional projection of the Sobol' sequence are nearly 
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identical, no improvement is seen. 

In the Brownian Bridge formulation, the Sobol' sequence with antithetic 
sampling is much better than either antithetic variables or Latin hypercube 
sampling. This suggests that it must be capturing some higher dimensional 
antisymmetric structure, most probably among the first few variables. 

The theory of scrambled (0, 360)-sequences in base 361 predicts that they 
will not be much better than Latin hypercube sampling until N = 361^ = 
130321 which is beyond the range we explore here. But by packing most of 
the structure into the first few dimensions of the integrand, we can consider 
methods in which scrambled nets are used on the most important dimensions 
and something else is used on the rest. Such methods have been proposed 
before: Spanier [19] describes a scheme in which quasi-Monte Carlo meth- 
ods are used on the first few dimensions of an integrand and simple Monte 
Carlo is used on the rest, and Owen [13] considered augmenting a random- 
ized orthogonal array with further dimensions taken from Latin hypercube 
samples. 

We considered using a scrambled (0, 32)-sequence in base 32 for the first 
32 dimensions and Latin hypercube samphng on the last 328 dimensions. 
For N = 32"^ one gets a scrambled (0, m, 32)-net in base 32 for the first 32 
dimensions and Latin hypercube sampUng for the rest. This will integrate 
the one dimensional part of the function and much of the m dimensional 
structure (superposition sense) of the first 32 dimensions, with the rest of the 
structure being integrated at the Monte Carlo rate. But (0, 32)-sequences 
can be stopped early or extended as necessary, whereas Latin hypercube 
sampling requires a prespecified number N of runs. As a compromise we ran 
a scrambled (0, 32)-sequence for the first 32 dimensions and took repeated 
independent Latin hypercube samples of size 1024 for the last 328 dimensions. 
Such a simulation can be conveniently stopped at any multiple of 1024 runs. 
The results are shown on Figure 4, labeled as RQR-BB for randomized quasi- 
random with Brownian bridge. All pairs of two variables among the first 32 
variables start to become balanced at sample size N = 32^ = 1024 and 
similarly, all triples of variables among the first 32 variables start to become 
balanced at sample size A'' = 32^ = 32768. This leads to results which are 
similar to the antithetic Sobol' with Brownian bridge results. However, the 
convergence rate for the scrambled net in the Brownian bridge representation 
appears to be larger, leading to greater accuracy at large N. 

We were able to achieve still better results than those shown here by 
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Figure 3: Error versus N (log base 10) for the Nearly Linear Problem, in the 
original representation. 

approximating the integrand by a quadratic function about the origin in 
the Gaussian coordinates and using this function as a control variate with 
antithetic random Monte Carlo. The quadratic terms were calculated by 
evaluating all 360 * 361/2 second derivatives of the integrand. 

6.2 Nonlinear Example 

We next consider the Nonlinear Example and display the same error curves. 
The mean value of PV for this problem is 130.712365. This value was ob- 
tained by averaging 200 independent Cranley-Patterson randomizations of 
2^^ fixed paths, generated by using the Sobol' sequence with the Brownian 
bridge ordering. These randomizations are due to Cranley and Patterson [2] 
and Tufhn [24] appears to be the first to realize their utility on nets. The 
certainty of this answer can be estimated by considering a six standard de- 
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Figure 4: Error versus N (log base 10) for the Nearly Linear Problem, in the 
Brownian bridge representation. 
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viation range, which was determined to be (130.712348,130.712382). The 
error curves in Figures 5 and 6 are based on using this value as the exsct 
solution. 

The mean length of a mortgage in this case is 76.5 months and the me- 
dian length is 58 months. The variance of PV in this value is 18.54 and 
the variance in the antithetic computation of this value is 1.127. Thus the 
function is about 97.7% antisymmetric. The variance under Latin hypercube 
sampling is 1.087 so that the function is only about 94.1% one dimensional. 
This may seem like a lot of one dimensional structure but compared with the 
previous example, the proportion of higher dimensional structure is greatly 
increased. 

As in the nearly linear example, antithetic sampling and Latin hypercube 
sampling in combination do not work better than separately. For this inte- 
grand as for that one, they each appear to remove the same source of error. 
In fact, for the nonlinear example, both give almost exactly the same errors. 

Figure 5 shows the error reference lines from Monte Carlo and antithetic 
Monte Carlo sampling. Both Latin hypercube sampling and the random- 
ized (0,360)-sequence in base 360, as well as their antithetic counterparts, 
give roughly the same accuracy as the antithetic Monte Carlo sampling. For 
simplicity, these error curves have therefore not been included in this graph. 
Superimposed are errors and lines for Sobol', and SoboP with antithetic sam- 
pling. In this case the Sobol' sequence is seen to catch up with the antithetic 
random sampling, but little is gained by combining antithetic sampling with 
the quasi-rajidom sequence. The Sobol' sequence outperforms Latin hyper- 
cube sampling on this, probably because the one dimensional parts of it are 
no longer so dominant. 

Figure 6 shows the errors from the Brownian bridge representation of the 
integrand. Reference lines are shown for Monte Carlo, antithetic sampling 
and Latin hypercube sampling. In this case also, Latin hypercube sampling 
does a bit better after the Brownian bridge transformation, suggesting that 
the function has become somewhat more one dimensional. As before, com- 
bining Latin hypercube sampling with antithetic sampling does not do much 
good. 

The Sobol' sequences perform especially well on the nonlinear function, 
in the BB representation with antithetic sampling. In terms of equation (3.4) 
this may be due to good equidistribution among the leading Sobol dimensions 
matched with their greater importance to the integrand. 
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7 Conclusions 



Our main conclusions are the following: 

• Quasi-Monte Carlo methods provide significant improvements in ac- 
curacy and computational speed for problems of small to moderate 
dimension. 

• While the effectiveness of quasi-Monte Carlo can be lost on problems of 
high dimension, this does not happen if the integrand is of low effective 
dimension in the superposition sense. 

• Some problems that have a large nominal dimension can be reformu- 
lated to have a moderate-sized effective dimension, so that the effec- 
tiveness of quasi-Monte Carlo is increased. 

• The Brownian bridge representation reduces the effective dimension for 
problems like the mortgage-backed security problem described here. 

Instead of straightforward use of high-dimensional quasi-random sequences, 
our recommendations are: 

• First analyze the problem, mathematically or numerically, to determine 
the most important input dimensions. 

• Where possible reformulate the problem to concentrate the variation 
in fewer dimensions. 

• When a small number of dominant dimensions can be identified or in- 
duced, apply quasi-random or randomized quasi-random (when sample 
based error estimates are desired) sequences to those dimensions. 

• For the remaining dimensions use pseudo-random or Latin hypercube 
sampling. 

• Consider applying classical variance reduction techniques, such as an- 
tithetic sampling, control variates, stratification and importance sam- 
pling in conjunction with the above. 
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We believe that high dimensional integration problems can range in diffi- 
culty from completely intractable to quite simple. In some cases it is possible 
to turn the former into the latter by carefully engineering the integrand. It 
is too early to say whether such manageable integrands are rare or dominant 
in financial applications. However, the results here indicate that, for the 
valuation of securities which depend on a single stochastic factor modeled 
as a Gaussian process with nonstochastic drift and volatility, the Brownian 
bridge representation may be extremely effective in reducing the dimension 
of the simulation. 
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